Largest Lyapunov Exponent for Many Particle Systems at Low Densities 
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The largest Lyapunov exponent A + for a dilute gas with short range interactions in equilibrium 
is studied by a mapping to a clock model, in which every particle carries a watch, with a discrete 
time that is advanced at collisions. This model has a propagating front solution with a speed that 
determines A + , for which we find a density dependence as predicted by Krylov, but with a larger 
prefactor. Simulations for the clock model and for hard sphere and hard disk systems confirm these 
results and are in excellent mutual agreement. They show a slow convergence of A + with increasing 
particle number, in good agreement with a prediction by Brunet and Derrida. 
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Recently, there has been great interest in the relation- 
ship between statistical mechanics and the theory of dy- 
namical systems 0-0] • Calculating dynamical properties 
such as Lyapunov exponents for statistical mechanical 
systems usually requires numerical simulations. For the 
Lorentz gas however, Dorfman, Van Beijeren and others 
have obtained analytical expressions for the Lyapunov 
spectrum and Kolmogorov-Sinai entropy at low densities, 
both in equilibrium and for the field-driven case ||,|J . 

In this paper we present an analytic calculation of the 
largest Lyapunov exponent in the low density limit for 
a gas at equilibrium consisting of particles with short 
range interactions. Our method is based on arguments 
from kinetic theory and similar in spirit to the method of 
Refs. We compare our results to those from com- 

puter simulations on hard disk and hard sphere systems 
and pay special attention to the dependence of the largest 
Lyapunov exponent on the total number of particles. 

We consider a gas consisting of N atoms of diameter a, 
defined as the (strictly finite) range of interaction, and 
mass m in d dimensions, in a volume V. The reduced 
density n is defined as Na d /V and will serve as a small 
parameter. To calculate the largest Lyapunov exponent 
we follow two nearby trajectories in phase space. For the 
first one, the reference trajectory, the positions and veloc- 
ities of the particles are denoted by (fj, Vi). In the second 
trajectory they are denoted by (fj + Sfi, Vi + SHi). The 
deviations (Sfi,Svi) will be taken to be infinitesimally 
small. For a chaotic system, they will grow exponen- 
tially with time at a rate equal to the largest Lyapunov 
exponent A + . Since the whole vector (Sri^Svi) in phase 
space grows exponentially, so will a generic projection, 
hence one has 



A + = lim — In 



E?=i\\sm\ 



Eli 11^(0)1 



(1) 



Therefore, in order to calculate A + one has to find out 
how 5vi{t) typically increases with time. We will first il- 
lustrate this on the somewhat simpler case of the random 



Lorentz gas, consisting of a single light particle moving 
among a random array of fixed scatterers interacting with 
the light particle through a spherically symmetric poten- 
tial. Between collisions the velocity deviation does not 
change and the position deviation changes according to 

5f%t) =5f f (t ) + (t-t )6v(t ). (2) 

In a collision the velocity changes from v to v given by 

v = v- 2(n • v)n = M. h v. (3) 

n denotes the unit vector in the direction from the cen- 
ter of the scatterer to the point of closest approach. The 
change of Sv in a collision is obtained from Eq. (|^) by 
expanding both v + Sv and n + Sh to linear order in the 
deviations. The difference in impact times for the two 
nearby trajectories leads to a shift in Sr. Since devia- 
tions follow linearized dynamics one always has 
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For hard sphere scatterers with radius a it turns out that, 
in any number of dimensions, A = B = M„, P = and 

Q = [u(n ■ v)]^ 1 [(n ■ v)l + nv\ ■ [(n ■ v)l — vn] , (5) 

with 1 the identity matrix. A derivation of these results 
in two dimensions can be found in ||. From the above 
equations we infer that at low density, just after the /c'th 
collision, with k very large, Sv and Sf will typically have 
increased to 



5v"{t k ) 



v (a/h) k 



S?(t k ) 



a(a/n) k . 
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with v the speed of the light particle, and a a constant 
of order unity. This follows from an inductive argu- 
ment: suppose Eq. (^) is valid after the fc'th collision, 
then according to Eqs. (g) and (9) one has Sr(tk+i) « 
Sf^tk) + t rn fSv'(tk) ~ a(a/n) k ^ , where we replaced 
tk+i — tk by its average value, the mean free time t m f. 
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In the last approximate equality we neglected Sf(tk) 
since it is one order of h smaller than t m fSv'(tk). Ac- 
cording to Eq. (H) and (0), after the {k + l)th collision 
5v'(tk+i) = 6v(tk) — 2QSr(tk+i) ~ v(a/h) h+1 , where we 
neglected Sv(tk) since it is one order of h smaller than the 
second term, and used that the typical size of the matrix 
elements of Q is (v/a), as is explicitly seen in Eq. (||). 
Now, because t « kt m f = k/v, with v the single particle 
collision frequency, it follows from Eqs. ([[]) and (|^) that 
the Lyapunov exponent is 

A + = —vh\n + v In a, 

with a to be determined by an averaging procedure over 
free flight times and collision dynamics. This estimate 
was already obtainted by Krylov Q. Notice that the 
value of a is not important for the dominant first term 
in A+. 

These considerations can be generalized to systems of 
identical moving particles by noting that in a collision, 
say between particles 1 and 2, Eqs. (|[^) still are appli- 
cable to the relative velocity, v = Vi — V2, the relative 
velocity deviation, Sv = Svi — Sv2 and the relative posi- 
tion deviation Sf — 5f\ — <5r*2 . In addition one needs the 
corresponding relations for the center of mass coordinates 
V = (vi + V2)/2 and R = (rj + f2)/2, which are 

V' = V ; 8V' = SV ; SR' = SR. (7) 

Assume now that the deviations for particles 1 and 2 
just after their last collisions before the present one were 
of the form (||) with exponents fe respectively fe, and 
v the mean relative velocity. By similar reasoning as for 
the Lorentz gas it follows that just before collision Sv and 
SV both are of order (a/n) max ( kl - k2 ' > whereas Sf and SR 
are of order (a/h) max< - kl ' k2 } +1 . As a consequence of (H) 
and (0) right after the collision Sf[ and Sv^ (i = 1, 2) will 
then also be of order (a/n)" la:E ^ fel ' fc2 ' +1 . So on average 
In \Sfi \ also increases by units of ln(a/n) at collisions, but 
in contrast to the Lorentz gas this increase may involve 
several of these units, in case the other particle involved 
in the collision has a higher fc-value. 

The values of In a in an actual realization of the dy- 
namics will fluctuate strongly from collision to collision. 
However, their distribution becomes independent of den- 
sity and increasingly narrow relative to ln(l/n) as density 
gets closer to zero. Therefore the essence of the dynamics 
determining the largest Lyapunov exponent is captured 
in the following simple clock model: Think of each parti- 
cle i as carrying a watch, whose clock value is fe. When 
two particles collide, they synchronize their watches to 
the larger of the two clock values, and advance them by 
one unit. The largest Lyapunov exponent will be de- 
termined by the speed w by which the watches run on 
average and will be of the form 

A + = w(— v\nn + ^lna). (8) 

The synchronization of the fc-values prohibits a direct 
identification with the number of collisions like we could 
do in the Lorentz gas. 



We will use a mean field approach to calculate the clock 
speed w. We denote the number of particles that have a 
given clock value k by Nk and assume that they are dis- 
tributed uniformly in V . In collisions involving particles 
with clock value k, Nk decreases. It is increased by two 
in collisions in which the largest incoming ki was k — 1. 
So the rate equations for the Nk become 

dN °° 

= - ~ 2i? (M) + 2 Z] R (k-l,t)- 

l= — oo l= — oo 

l^k 

R(k,i) are the rates by which collisions between k and I 
take place. We use a StoBzahlansatz: the rate of collisions 
between particles with clock values k and I is propor- 
tional to NkNi/N 2 . Since all rates are also proportional 
to v 1 we will express time in units of the mean free time: 
t = vt. We use the fractions fk = Nk/N to eliminate 
the TV dependence: 

-A = ~fk+ 2/fe-l ^2 fl + fk-V 



For the cumulatives Ck — J2i=-oo fi this reduces to 

^+C k = C 2 k _ v (9) 

The solution is given by the recursion relation 

C k {r) = e- T C k (0) + f e T '- T Cl_ x (T')dT' . 
Jo 

If Ck is zero at r = it remains zero. Thus the start- 
ing point of this recursion is the smallest k for which 
Cfe(O) ^ 0. Inductively we see that all C k are polynomials 
in e~ r , of which the degree grows exponentially with k. 
We calculated these polynomials with initial conditions 
corresponding to fk{r = 0) = Ski- The exponentially 
growing degree of the polynomials enables only a limited 
number of Ck to be computed, even on a computer. The 
results up to k = 30 at several time values are shown in 
Fig. [I]. The initial distribution broadens and moves to 
the right. We expect the distribution to asymptotically 
become a front propagating at a constant speed w. Then 
we have 

N oo 

£n^(i)ii 2 = £ am , 2 e - 2fel ^) 

i— 1 k— — oo 

oo 

k— — oo 
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FIG. 1. The cumulative distribution of clock values k at 
times t — 0, t — 1, r = 2 . . .r = 10 (from left to right). 

This result should be proportional to e 2tA+ so one indeed 
recovers Eq. (||). It agrees with Krylov's conjecture 
except for the appearance of w. The result of Stoddard 
and Ford [M is of the same form when expanded in n. 

The speed w should be independent of the initial condi- 
tions. Physical initial conditions will have finite support, 
because there is always a particle with smallest k and a 
particle with highest k. Then Cfc(O) is bounded by two 
step functions. Using Eq. @ one can show that they will 
remain bounded by the solutions corresponding to these 
two initial conditions. If these tend to some uniformly 
moving solution with speed w, so will the real system. 
Thus w is unique for this set of initial conditions. 

We put the propagating front Ansatz Cfe(r) = F(k — 
wt) into Eq. (^|) to obtain a differential-difference equa- 
tion M for the shape of the cumulatives: 



dF 

■ w—(x) + Fix) = F 2 (x - 1). 
ax 



(10) 



F has to be monotonically increasing, tending to as 
x — > — oo and to 1 as x — > oo. This means that F = has 
to be unstable and F = 1 has to be stable. It is easy to see 
that these are fixed points of Eq. (|l^) . Their stability is 
determined by linearized equations. The behavior around 
a fixed point is always exponential: F(x) — J2jPj eSjX > 
in which the Sj are roots of the so-called characteris- 
tic equation and pj are polynomials in x of degree less 
than the multiplicity of root Sj || . For an unstable fixed 
point some of the Sj have positive real parts. Around 
F = 0, this is true if w > 0. For a stable fixed point, 
the term with least negative Sj, let's call this —7, will 
dominate the large x behavior. If 7 were complex, we 
would see oscillatory behavior, violating monotonicity, 
so 7 has to be real. Inserting the asymptotic behav- 
ior F(x) = 1 — exp(— jx) into Eq. ( |l0| ) and neglect- 
ing quadratic terms, produces the characteristic equation 
71/7+ 1 — 2e 7 = 0. This gives a relation between w and 7: 
w (l) = (2e 7 — 1)/7- It turns out that there is a minimum 
w for positive real values of 7 which can be expressed in 
terms of Lambert's W function: 



w = -l/W(-l/2e) = 4.31107.. 

This value is in accordance with estimates from Fig. [I]. 
Solutions with initial conditions with finite support select 
this minimum speed. The same kind of velocity selection 
occurs in other systems, for a number of which it has 
been proved ||. 

We compared the result w = 4.311.. with those from 
simulations done by Dellago, Posch and Hoover ||, with 
64 hard disks. They made a fit of the largest Lyapunov 
exponent to a nhxin/b) indicating a value of w rj 3.3. 
The difference to our value turns out to be due to large 
finite A effects. First we'll show this numerically. We 
take A watches and give them some initial k values. In 
each time step we pick two watches at random and ad- 
vance their fc-values according to the rules of the clock 
model. We compute the average growth of k per watch 
per time step, and find wn- We did this for numbers of 
watches ranging from 4 to 2 19 . In Fig. § the results were 
fitted to an algebraic curve 



w N 



4.311 -AN~ 



(11) 



A good fit, except for the smallest values of A, was ob- 
tained by choosing B = 0.277 and A = 3.466. No good 
fit could be obtained for B — 1 or on replacing the al- 
gebraic A-dependence by an exponential one. The value 
B = 0.277 is in reasonable agreement with the exponent 
of —1/3 obtained by Dellago and Posch pO] . In order to 
come to a better comparison between clock model predic- 
tions for wn and simulation results on actual dilute gas 
models we performed a number of new simulations using 
the same methods as in Refs. [| 10 , both on hard disk 
and hard sphere systems for different particle numbers 
at a number of low densities, n = 10~ 5 , 10~ 4 , 10~ 3 and 
10 -2 . For each N the results for the largest Lyapunov 
exponent were fitted to — wuhxin/a). The results for w 
are also plotted in Fig. ^. One sees that the results of 
the clock model are in excellent agreement with those of 
the simulations. 

Stoddard and Ford used crude arguments to get 
wjv = In A. This relation is also plotted in Fig. |2[ One 
sees that this only gives a good fit for very small N. 
In simulations of 100 particles with a cut-off Lennard- 
Jones interparticle potential, Stoddard and Ford found 
agreement with their predicted value, which lies some- 
what above our asymptotic value of 4.311 and much more 
above the simulation results for 100 particles, both in the 
clock model and for hard spheres and disks. Stoddard 
and Ford acknowledged that the error increases as n gets 
smaller and say that their simulation results for low n 
should not be expected to fit their theory. But the n at 
which the error becomes too big, is not sharply defined. 
If one takes it low enough, the data support their predic- 
tion, but if one takes it a little higher the results come 
more in line with those from the clock model simulations. 

In a recent paper pd| , Brunet and Derrida present a 
way to compute the A dependence of the velocity in a 
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similar model by treating it as a discretization effect. In 
our case, there are at least two particles with highest k 
in any realization. This means we have to take e = 2/N 
in equation (7) in flff . Inserting our expression for w(-/), 
we find 



W N 



l)ir 2 



2W{N/2)' 



(12) 



We plotted this prediction also in Fig. ^. The agreement 
is good for N > 100. 

In the work by Searles et al |L2| a weak but persistent 
increase in A + with N was interpreted as a sign of a log- 
arithmic divergence. It was argued that the data were 
not consistent with a 1/A-approach to a constant value 
and a plot of A + versus lnA^ looks quite linear over an 
appreciable range. However, Dellago and Posch jl(| in 
their simulations on dense hard sphere systems did not 
observe such a divergence and in fact it looks like the re- 
sults of Searles et al are entirely consistent with the type 
of behavior predicted by Brunet and Derrida. The mean 
field analysis given here is not decisive though, since it 
completely ignores all effects of local density and temper- 
ature fluctuations. 



ranged, i.e. it strictly vanishes beyond its range a, or per- 
haps may be allowed to vanish exponentially. This let- 
ter shows that the calculation of dynamical properties of 
many particle systems is feasible and that the calculation 
of the largest Lyapunov exponent in dilute gases requires 
the solution of a nonlinear front propagation equation. 
The method will be extended in future work to get the 
O(h) term of the Lyapunov exponent. This term will de- 
pend on the details of the interaction at a collision, and 
is of considerable physical interest. 
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FIG. 2. Number dependent coefficient w. The circles are 
clock model results. The bars are new molecular simulation 
results: wide error bars for hard spheres, narrow error bars 
for hard disks. The dashed line is Stoddard and Ford's lniV 
prediction. The dashed-dotted curve is a fit of the mean field 
results to the algebraic expression (|ll|). The thick line gives 
the analytic result for N — > oo. The solid curve is the predic- 
tion of Eq. (||). 

We conclude by stressing that the first term of the 
density expansion of the largest Lyapunov exponent of 
a dilute gas that was calculated in this letter, is univer- 
sal for systems where the interaction is sufficiently short 
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